Output contents

Notes about this data

Set Gobal Chunk Options & Load libraries

knitr::opts_chunk$set(
    echo = T,
    warning = FALSE,
    cache = FALSE,
    collapse = T
)
getwd()
library(tidyverse) #WHY WAS THIS NOT LOADED PREVIOUSLY??!!!
#library(nlme) #mixed modeling - delete if not used
library(plyr)#revalue() factors and reorder them! dplyr::recode not working
library(dplyr)  #summarise_at()
library(ggplot2) #graphing
#library(lattice) #delete if not used
#library(gtable)
library(grid) #tables?
library(gridExtra)#masks combine() from dplyr
library(pander) #table output
library(Publish) #delete if not used
library(ggpubr) #combine plots
library(kableExtra)
library(devtools) #
library(broom) #
library(purrr)#
#Modeling libraries
library(car)#sophisticated anova() modeling
#library(rstatix)#anova_test() - masks mutate() and desc() from plyr package!
library(sjstats)#anova_stats()
library(pwr)#unknown function but the sjstats package must use it to calculate Cohen's D!
#library(MDmisc)#download from GitHub - 

Import Data & Set Factor Levels

#rm(Recovery)
D<-read.csv("T0_Stress_R.csv")

# Set factors (no longer automatic)
D$Moss<-as.factor(D$Moss)
D$Microsite<-as.factor(D$Microsite)
D$Habitat<-as.factor(D$Habitat)

#Remove other species & drop PAR rejected for tissue sampling
D<-droplevels(D[!c(D$Microsite=="PJ-N-91" | D$Microsite =="PJ-S-67a" | 
                    D$Microsite=="BB-N-35" | D$Microsite=="BB-S-53"),])
length(levels(D$Microsite))#Down to 94 units

# Order factor levels
D$Topography<-factor(D$Topography, levels=c("S", "F", "N"))
D$Veg.Type<-factor(D$Veg.Type, levels=c("CB", "BB", "PJ"))
D$Veg.Type<-plyr::revalue(D$Veg.Type, c("CB" = "Low-scrubland", "BB" = "Mid-shrubland", "PJ" ="High-woodland"))

D$Site<- factor(D$Site, levels=c("1","2","3"))
#levels(D$Site)<- list("1" = "Low-scrubland", "2"= "Mid-shrubland", "3" = "High-woodland")
## [1] 94

Script error correction (PhiPSII, qP)

#Mean PhiPSII (Phi-PSII) of the Assays (factor levels)
D %>% group_by(., Assay) %>% summarise_at(vars(PhiPSII), mean)
#Mean difference (Assay 1 - the others)
with(D, mean(PhiPSII[Assay!=1]))#0.2464286
## [1] 0.2464286
with(D, mean(PhiPSII[Assay!=1])) - with(D, mean(PhiPSII[Assay==1]))#-0.1834881 
## [1] -0.1834881
#D<-D %>% mutate(PhiPSII = ifelse(Assay ==1, PhiPSII-0.17, PhiPSII))# Correction = -0.183
#D$Moss[D$PhiPSII< 0]#Chosen correction where no values become negative in Assay 1.
#D<-D %>% mutate(qPhoto = ifelse(Assay ==1, qP-0.22, qP))
#D$Moss[D$qPhoto<0]#Worked all are positive

#Correlation of Fv/Fm for Assay 2 - 4
#Assay 2 = 0.9056776, Assay 3 = 0.9093973, Assay 4 = 0.9148743, Assay 1 = 0.9627869
D %>% group_by(., Assay) %>% do(as.data.frame(cor(.$Fv.Fm, .$PhiPSII), data=.))

#Correlation of the pooled dataset, thus, improves after correction, but too much?
with(D, cor(PhiPSII, Fv.Fm))#0.74
## [1] 0.7444226
#with(D, cor(PhiPSII, Fv.Fm))#0.909 

Delete Assay 1

#See change after correction medians align!
par(mar = c(4, 4, .1, .1))
with(D, boxplot(PhiPSII~Assay))
with(D[D$Assay != 1,], boxplot(PhiPSII~Assay))#Dropped Assay 1
# Density plots;transparent
ggplot(D[D$Assay != 1,], aes(x=PhiPSII, fill=Veg.Type)) + geom_density(alpha=.6)#

#————————- # SUMMARY STATS ## Function * New way to enter functions after 2019 packages update:

summary(D)# means, min, max, etc for all variables
##       Time        Moss     Moss.number      Moss.sub            Microsite 
##  Min.   :0   10     : 1   Min.   : 1.00   Length:94          BB-F-37 : 1  
##  1st Qu.:0   11     : 1   1st Qu.:22.25   Class :character   BB-F-38a: 1  
##  Median :0   12     : 1   Median :47.50   Mode  :character   BB-F-39a: 1  
##  Mean   :0   13a    : 1   Mean   :47.12                      BB-F-40a: 1  
##  3rd Qu.:0   14a    : 1   3rd Qu.:71.75                      BB-F-41a: 1  
##  Max.   :0   15a    : 1   Max.   :96.00                      BB-F-42a: 1  
##              (Other):88                                      (Other) :88  
##  Old.Habitat              Habitat       Index       Site            Veg.Type 
##  Length:94          Canopy    :48   Min.   : 3.00   1:26   Low-scrubland:26  
##  Class :character   Dripline  :28   1st Qu.: 8.00   2:34   Mid-shrubland:34  
##  Mode  :character   Interspace:18   Median :14.00   3:34   High-woodland:34  
##                                     Mean   :14.44                            
##                                     3rd Qu.:20.00                            
##                                     Max.   :26.00                            
##                                                                              
##  Topography  Shade.Index         PDIR            Slope            Aspect      
##  S:23       Min.   : 6.00   Min.   :0.7717   Min.   : 0.000   Min.   :  0.00  
##  F:37       1st Qu.:12.00   1st Qu.:0.9142   1st Qu.: 2.000   1st Qu.: 90.75  
##  N:34       Median :20.00   Median :0.9846   Median : 7.000   Median :200.00  
##             Mean   :17.87   Mean   :0.9486   Mean   : 7.468   Mean   :194.88  
##             3rd Qu.:23.00   3rd Qu.:0.9997   3rd Qu.:12.000   3rd Qu.:270.00  
##             Max.   :27.00   Max.   :1.0496   Max.   :18.000   Max.   :357.00  
##                                                               NA's   :4       
##    VegPercent        SolarVeg          Camera           Pair     
##  Min.   :  0.00   Min.   :  0.00   Min.   : 1.00   Min.   : 1.0  
##  1st Qu.: 22.22   1st Qu.: 29.17   1st Qu.: 6.25   1st Qu.: 9.5  
##  Median : 38.89   Median : 43.75   Median :11.50   Median :19.0  
##  Mean   : 40.66   Mean   : 45.26   Mean   :11.73   Mean   :18.9  
##  3rd Qu.: 60.41   3rd Qu.: 62.50   3rd Qu.:16.75   3rd Qu.:28.5  
##  Max.   :100.00   Max.   :100.00   Max.   :23.00   Max.   :37.0  
##                                    NA's   :72      NA's   :55    
##  Notes.Urgent        Assay.Set             Assay            Clip      
##  Length:94          Length:94          Min.   :1.000   Min.   : 1.00  
##  Class :character   Class :character   1st Qu.:1.250   1st Qu.: 6.25  
##  Mode  :character   Mode  :character   Median :2.000   Median :12.00  
##                                        Mean   :2.468   Mean   :12.32  
##                                        3rd Qu.:3.000   3rd Qu.:18.00  
##                                        Max.   :4.000   Max.   :24.00  
##                                                                       
##        Fo               Fm              Fv             Fv.Fm       
##  Min.   : 30.00   Min.   : 52.0   Min.   : 12.00   Min.   :0.1300  
##  1st Qu.: 79.25   1st Qu.:129.5   1st Qu.: 45.25   1st Qu.:0.3515  
##  Median : 96.00   Median :183.5   Median : 93.50   Median :0.4770  
##  Mean   : 99.02   Mean   :195.0   Mean   : 95.97   Mean   :0.4546  
##  3rd Qu.:118.50   3rd Qu.:248.8   3rd Qu.:132.00   3rd Qu.:0.5697  
##  Max.   :174.00   Max.   :419.0   Max.   :265.00   Max.   :0.7370  
##                                                                    
##        Fs              Fm.L            Fo.L             Fv.L       
##  Min.   : 29.00   Min.   : 46.0   Min.   : 17.00   Min.   : 15.00  
##  1st Qu.: 76.25   1st Qu.:101.0   1st Qu.: 52.50   1st Qu.: 42.25  
##  Median : 92.50   Median :137.5   Median : 71.50   Median : 61.00  
##  Mean   : 96.63   Mean   :141.8   Mean   : 70.89   Mean   : 70.86  
##  3rd Qu.:118.00   3rd Qu.:175.5   3rd Qu.: 89.25   3rd Qu.: 95.00  
##  Max.   :190.00   Max.   :319.0   Max.   :149.00   Max.   :211.00  
##                                                                    
##     Fv.Fm.L          PhiPSII             qP              qNP         
##  Min.   :0.2230   Min.   :0.0430   Min.   :0.1330   Min.   :-0.8570  
##  1st Qu.:0.4050   1st Qu.:0.1900   1st Qu.:0.4903   1st Qu.: 0.1765  
##  Median :0.4905   Median :0.2685   Median :0.5820   Median : 0.4100  
##  Mean   :0.4881   Mean   :0.2933   Mean   :0.5826   Mean   : 0.3651  
##  3rd Qu.:0.5780   3rd Qu.:0.3635   3rd Qu.:0.6707   3rd Qu.: 0.5870  
##  Max.   :0.7360   Max.   :0.6930   Max.   :0.9460   Max.   : 0.7970  
##                                                                      
##       NPQ             Shrubindex        Topoindex        Siteindex       
##  Min.   :-0.19100   Min.   :-1.0000   Min.   :-1.000   Min.   :-1.00000  
##  1st Qu.: 0.08575   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.:-1.00000  
##  Median : 0.32000   Median : 1.0000   Median : 0.000   Median : 0.00000  
##  Mean   : 0.43293   Mean   : 0.3191   Mean   : 0.117   Mean   : 0.08511  
##  3rd Qu.: 0.68800   3rd Qu.: 1.0000   3rd Qu.: 1.000   3rd Qu.: 1.00000  
##  Max.   : 2.54200   Max.   : 1.0000   Max.   : 1.000   Max.   : 1.00000  
##                                                                          
##   Multibuffer     
##  Min.   :-2.0000  
##  1st Qu.: 0.0000  
##  Median : 1.0000  
##  Mean   : 0.5213  
##  3rd Qu.: 1.0000  
##  Max.   : 3.0000  
## 

#Stats needed for plotting
SE<- function(x) round((sqrt(var(x)/length(x))),3)
Myfunctions = list (
  Mean = function(x) round(mean(x),3),
  Max= function(x)round(max(x),3),
  Min= function(x)round(min(x),3),
  Range= function(x) round((max(x)-min(x)),3),
  SE= function(x) round((sqrt(var(x)/length(x))),3),
  Upper= function(x) round(mean(x)+SE(x),3),
  Lower= function(x) round(mean(x)-SE(x),3),
  Var = function(x) round(var(x),3),
  SD = function(x) round(sd(x),3)
)

Tables

###Stress (Fv/Fm) & Performance (PhiPSII)

#Fv/Fm Tables
#Site table
T.Site.Fvfm<-D %>% 
  group_by(Site) %>% 
  summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Site.Fvfm, caption = "Fv/Fm Stats by Site", format="pandoc", padding=2)
Fv/Fm Stats by Site
Site Mean Max Min Range SE Upper Lower Var SD
1 0.440 0.602 0.183 0.419 0.021 0.461 0.419 0.011 0.105
2 0.429 0.737 0.130 0.607 0.027 0.456 0.402 0.026 0.160
3 0.491 0.686 0.144 0.542 0.028 0.519 0.463 0.026 0.163
#Topography tables
T.Topo.Fvfm<-D %>% 
  group_by(Topography) %>% 
  summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Topo.Fvfm, caption = "Fv/Fm Stats by Topography", format="pandoc", padding=2)
Fv/Fm Stats by Topography
Topography Mean Max Min Range SE Upper Lower Var SD
S 0.432 0.656 0.20 0.456 0.026 0.458 0.406 0.016 0.125
F 0.391 0.647 0.13 0.517 0.027 0.418 0.364 0.026 0.162
N 0.539 0.737 0.25 0.487 0.018 0.557 0.521 0.012 0.107
#Microhabitat tables
T.Hab.Fvfm<-D %>% 
  group_by(Habitat) %>% 
  summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Hab.Fvfm, caption = "Fv/Fm Stats by Habitat", format="pandoc", padding=2)
Fv/Fm Stats by Habitat
Habitat Mean Max Min Range SE Upper Lower Var SD
Canopy 0.516 0.737 0.174 0.563 0.020 0.536 0.496 0.018 0.136
Dripline 0.421 0.647 0.130 0.517 0.028 0.449 0.393 0.021 0.146
Interspace 0.342 0.508 0.160 0.348 0.025 0.367 0.317 0.011 0.105

#Phi-PSII Tables
##Site Table
T.Site.Phi<-D %>% 
  group_by(Site) %>% 
  summarise_at(vars(PhiPSII), .funs=Myfunctions)
kable(T.Site.Phi, caption = "PHi-PSII Stats by Site", format="pandoc", padding=2)
PHi-PSII Stats by Site
Site Mean Max Min Range SE Upper Lower Var SD
1 0.256 0.520 0.150 0.37 0.019 0.275 0.237 0.010 0.098
2 0.281 0.693 0.043 0.65 0.024 0.305 0.257 0.020 0.141
3 0.334 0.622 0.072 0.55 0.026 0.360 0.308 0.023 0.151
#Topo table
T.Topo.Phi<-D %>% 
  group_by(Topography) %>% 
  summarise_at(vars(PhiPSII), .funs=Myfunctions)
kable(T.Topo.Phi, caption = "PHi-PSII Stats by Topography", format="pandoc", padding=2)
PHi-PSII Stats by Topography
Topography Mean Max Min Range SE Upper Lower Var SD
S 0.301 0.585 0.098 0.487 0.026 0.327 0.275 0.016 0.127
F 0.238 0.621 0.043 0.578 0.020 0.258 0.218 0.015 0.122
N 0.348 0.693 0.072 0.621 0.024 0.372 0.324 0.019 0.139
#Microhabitat table
T.Hab.Phi<-D %>% 
  group_by(Habitat) %>% 
  summarise_at(vars(PhiPSII), .funs=Myfunctions)
kable(T.Hab.Fvfm, caption = "Phi-PSII Stats by Habitat", format="pandoc", padding=2)
Phi-PSII Stats by Habitat
Habitat Mean Max Min Range SE Upper Lower Var SD
Canopy 0.516 0.737 0.174 0.563 0.020 0.536 0.496 0.018 0.136
Dripline 0.421 0.647 0.130 0.517 0.028 0.449 0.393 0.021 0.146
Interspace 0.342 0.508 0.160 0.348 0.025 0.367 0.317 0.011 0.105
#————- ——— –
# BOXPLOTS & M EANS

Life zone site boxplots

  • Outlier (Phi-PSII) is Moss 34 at Site 2. How could this reading be SO high? Was Fv/Fm high also? Yes: 0.737.
  • Question: what dimensions to save for publication? (fig.width=4, fig.height = 1.7)
#Fv/Fm life zone boxplots
PlotA<-ggplot(D, aes(x=(Site), y=Fv.Fm))+
  geom_boxplot(outlier.shape=NA, aes(x=factor(Site), fill= (Site)))+
  geom_point(aes(fill=(Site)), position = position_dodge2(width=0.3), size=1, 
             color="black")+
  #stat_summary(fun.y=mean, geom="point", shape=20, size=3, color="white")+
  geom_point(data = T.Site.Fvfm, aes(x=Site, y=Mean), color="white")+
  geom_errorbar(data=T.Site.Fvfm, aes(x=Site, y= Mean, ymin=Lower, ymax=Upper), width=.2, size=0.7, color="white") +
  theme_linedraw()+
  labs(y='\u2190 Stress (Fv/Fm)', x="Elevation life-zone (site)")+
  #scale_x_discrete(labels=c("1" = "Low-scrubland", "2" = "Mid-shrubland", "3"= "High-woodland"))+
  scale_x_discrete(labels=c("1" = "Low-site", "2" = "Mid-site", "3"= "High-site"))+
  scale_y_continuous(breaks = seq(0, 0.75, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75))+
  theme(legend.position = "none", text=element_text(size=15), axis.title.x = 
          element_blank(), panel.grid.minor = element_blank())+
  scale_fill_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"))# panel.grid.major = element_blank(),
#PlotA

#Phi-PSII life zone boxplots
PlotB<-ggplot(D, aes(x=(Site), y=PhiPSII))+
  geom_boxplot(outlier.shape=NA, aes(x=factor(Site), fill= (Site)))+
  geom_point(aes(fill=(Site)), position = position_dodge2( width=0.3), size=1, 
             color="black")+
  #stat_summary(fun.y=mean, geom="point", shape=20, size=3, color="white")+
  geom_point(data = T.Site.Phi, aes(x=Site, y=Mean), color="white")+
  geom_errorbar(data=T.Site.Phi, aes(x=Site, y= Mean, ymin=Lower, ymax=Upper), width=.2, size=0.7, color="white") +
  theme_linedraw()+
  labs(y='Performance (\u03A6PSII)', x="Elevation life-zone (site)")+
  scale_x_discrete(labels=c("1" = "Low-site", "2" = "Mid-site", "3"= "High-site"))+
  scale_y_continuous(breaks = seq(0, 0.75, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75))+
  theme(legend.position = "none", text=element_text(size=15), axis.title.x = 
          element_blank(), panel.grid.minor = element_blank())+
  scale_fill_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"))# panel.grid.major = element_blank(),
#PlotB

#Plot together
PlotsAB <-ggarrange(PlotA, PlotB, labels=c("a", "b"), widths=1:1, ncol=2, nrow=1, font.label= list(size=15)) #Add labels, change text sizes
PlotsAB<- annotate_figure(PlotsAB, top = text_grob("Summer stress assay vs. elevation-life zone sites", size = 15, rot=0, just= "center")) #Add titles
PlotsAB

Topography Boxplots

#Fv/Fm topo boxplots
PlotC<-ggplot(D)+
    aes(x=Topography, y = Fv.Fm, fill = Topography) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Topo.Fvfm, aes(x=Topography, y= Mean, ymin=Lower, ymax=Upper), 
                width=.2, size=0.7, color="white") +
    scale_y_continuous(name = "\u2190 Stress signal (Fv/Fm)",breaks = seq(0, 0.7, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75)) +
    scale_x_discrete(name = "Topography (pooled life zones)", 
                     labels=c("S-facing", "Flat", "N-facing")) +
   #ggtitle("Summer stress assay (Series 1.1)") +
      theme_linedraw() +
      theme(text = element_text(size=15), legend.position = "none") +
         scale_fill_manual(values = c("#fe9929", "#666666","darkorchid"))+
       theme(panel.background = element_rect(fill = "white"), axis.title.x = element_blank(), panel.grid.minor = element_blank())
#PlotC

#Phi topo boxplots
PlotD<-ggplot(D)+
    aes(x=Topography, y = PhiPSII, fill = Topography) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Topo.Phi, aes(x=Topography, y= Mean, ymin=Lower, ymax=Upper), 
                width=.2, size=0.7, color="white") +
    scale_y_continuous(name = "Performance (\u03A6PSII)",breaks = seq(0, 0.7, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75)) +
    scale_x_discrete(name = "Topography (pooled life zones)", 
                     labels=c("S-facing", "Flat", "N-facing")) +
   #ggtitle("Summer stress assay (Series 1.1)") +
      theme_linedraw() +
      theme(text = element_text(size=15), legend.position = "none") +
         scale_fill_manual(values = c("#fe9929", "#666666","darkorchid"))+
       theme(panel.background = element_rect(fill = "white"), axis.title.x = element_blank(), panel.grid.minor = element_blank())
#PlotD

#Plot together
PlotsCD <-ggarrange(PlotC, PlotD, labels=c("c", "d"), widths=1:1, ncol=2, nrow=1, font.label= list(size=15)) #Add labels, change text sizes
PlotsCD<- annotate_figure(PlotsCD, top = text_grob("Summer stress assay vs. Topography zones", size = 15, rot=0, just= "center")) #Add titles
PlotsCD

Microhabitat boxplots

#Fv/fm
PlotE<-ggplot(D)+
    aes(x=Habitat, y =Fv.Fm, fill = Habitat) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Hab.Fvfm, aes(x=Habitat, y= Mean, ymin=Lower, ymax=Upper), 
                width=.15, size=0.7, color= "white") +
    scale_y_continuous(name = "\u2190 Stress signal (Fv/Fm)",breaks = seq(0, 0.75, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75)) +
    scale_x_discrete(name = "Microhabitat type") +
  # ggtitle("Micro-habitat buffering") +
      theme_linedraw() +
      theme(text = element_text(size=15), legend.position = "none") +
         scale_fill_manual(values = c("darkgreen", "yellowgreen", "tan"))+
       theme(panel.background = element_rect(fill = "white"), panel.grid.minor = element_blank(),axis.title.x = element_blank())
#PlotE 

#Phi
PlotF<-ggplot(D)+
    aes(x=Habitat, y =PhiPSII, fill = Habitat) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.3), size=1, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Hab.Phi, aes(x=Habitat, y= Mean, ymin=Lower, ymax=Upper), 
                width=.15, size=0.7, color= "white") +
    scale_y_continuous(name = "Performance (\u03A6PSII)", breaks = seq(0, 0.75, 0.1), 
                       expand=c(0,0), limits=c(-0.01,0.75)) +
    scale_x_discrete(name = "Microhabitat type") +
  # ggtitle("Micro-habitat buffering") +
      theme_linedraw() +
      theme(text = element_text(size=15), legend.position = "none") +
         scale_fill_manual(values = c("darkgreen", "yellowgreen", "tan"))+
       theme(panel.background = element_rect(fill = "white"), panel.grid.minor = element_blank(),axis.title.x = element_blank())
#PlotF 

#Plot together
PlotsEF <-ggarrange(PlotE, PlotF, labels=c("e", "f"), widths=1:1, ncol=2, nrow=1, font.label= list(size=15)) #Add labels, change text sizes
PlotsEF<- annotate_figure(PlotsEF, top = text_grob("Summer stress assay vs. Habitat types", size = 15, rot=0, just= "center")) #Add titles
PlotsEF

Final boxplots: Stress & performance vs 3 scales of buffering

#Plot A-F together
PlotsA_F <-ggarrange(PlotA, PlotB, PlotC, PlotD, PlotE, PlotF, labels=c("a","b","c","d","e", "f"), widths=1:1, ncol=2, nrow=3, font.label= list(size=20)) #Add letters & up the size
PlotsA_F<- annotate_figure(PlotsA_F, top = text_grob("Summer stress & performance vs. 3 scales of habitat buffering", size = 16, rot=0, just= "center")) #Add bottom/top titles
PlotsA_F

#——————— # VARIANCE TESTS

  • Goal: Test null of equal variance across categorical levels. Perform batch analysis.
  • Test: Fligner-Killeen homogeneity of variance test: stats::fligner.test(). For when groups lack normal distributions.
  • Results: Only Fv/Fm ~ Topography is heteroskedastic (for the photosynthesis Anova’s; ignoring vegetation shade metrics for now).
#batch variance tests for life zone sites (Veg.Type)
V.tests<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(~ tidy(fligner.test(. ~ D$Veg.Type)),.id="Variable")#put results into a tidy df
kable(V.tests, format = "pandoc", caption = "Site Variance: Fligner-Killeen Equal Variance by Site: T0 summer stress (n=94)", digits= c(0, 3, 4))#plot pretty & round output

#Topography
V.tests<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(~ tidy(fligner.test(. ~ D$Topography)),.id="Variable")#put results into a tidy df
V.tests<-V.tests %>% select(-parameter, -method) #exclude output parameter & method

#plot pretty
kable(V.tests, format = "pandoc", caption = "Topography Variance: Fligner-Killeen Equal Variance by Topozone: T0 summer stress (n=94)", digits= c(0, 3, 4))#round output

#Microhabitat
V.tests<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(~ tidy(fligner.test(. ~ D$Habitat)),.id="Variable")#put results into a tidy df
V.tests<-V.tests %>% select(-parameter, -method) #exclude output parameter & method

#plot pretty
kable(V.tests, format = "pandoc", caption = "Habitat Variance: Fligner-Killeen Equal Variance by Microhabitat type: T0 summer stress (n=94)", digits= c(0, 3, 4))#round output
Site Variance: Fligner-Killeen Equal Variance by Site: T0 summer stress (n=94)
Variable statistic p.value parameter method
Fv.Fm 5.289 0.0711 2 Fligner-Killeen test of homogeneity of variances
PhiPSII 5.627 0.0600 2 Fligner-Killeen test of homogeneity of variances
Topography Variance: Fligner-Killeen Equal Variance by Topozone: T0 summer stress (n=94)
Variable statistic p.value
Fv.Fm 8.398 0.0150
PhiPSII 1.271 0.5297
Habitat Variance: Fligner-Killeen Equal Variance by Microhabitat type: T0 summer stress (n=94)
Variable statistic p.value
Fv.Fm 1.498 0.4727
PhiPSII 4.100 0.1287

#———————— # ONE-WAY ANOVA’s

  • Goal: use a family of univariate 1-way ANOVA’s to test for mean differences in photosynthethetic summer stress (Fv/Fm & Phi-PSII) across categorical predictors: elevation-life zone sites, topography zones (aka topozones), & microhabitat types.
  • Unbalanced Design Challenge: missing factor levels (South-slope at Site 1 and interspaces at Site 3) prevent a global 3-way ANOVA. Also, all levels have different sample sizes for all three predictors (another type of unbalanced design).
  • Heteroskedasticity Challenge: only Fv/Fm by Topozones is heteroskedastic. For this test I will use Welch’s ANOVA if assumption of normality is met. Otherwise, classic ANOVA is used (whether or not normality is perfect).
  • P-values: I explicitly avoid adjusting pvalues for this family of tests because this is exploratory analysis using messy ecological data. We want to report all possible patterns of significance!

Classic ANOVAs

#1-way classic ANOVAs - use var.equal = TRUE
W.Aovs<- D %>% select(Fv.Fm, PhiPSII) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Veg.Type, var.equal= TRUE, data= D)), 
                .id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, denom.df
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Veg.Type, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Veg.Type, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Classic ANOVAs vs Site - assumes equal variance", digits = c(1, 1, 3, 3, 4, 5, 3), table.attr = "style = \"color: black;\"") %>% kable_styling()
Classic ANOVAs vs Site - assumes equal variance
Variable num.df den.df statistic p.value Cohens_F Rsq
Fv.Fm 2 91 1.614 0.2048 0.188 0.034
PhiPSII 2 91 2.702 0.0724 0.244 0.056

## Multiple Comparisons
Mod<- aov(PhiPSII ~ Veg.Type, data = D, var.equal= TRUE)
TukeyHSD(Mod)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = PhiPSII ~ Veg.Type, data = D, var.equal = TRUE)
## 
## $Veg.Type
##                                   diff          lwr       upr     p adj
## Mid-shrubland-Low-scrubland 0.02457240 -0.058961924 0.1081067 0.7635722
## High-woodland-Low-scrubland 0.07795475 -0.005579571 0.1614891 0.0725895
## High-woodland-Mid-shrubland 0.05338235 -0.024383869 0.1311486 0.2360452

Welch’s ANOVAs: Topography

  • use stats::oneway.test() with argument var.equal = FALSE, for Welch Correction, an approximate method of Welch (1951), which generalizes the 2-sample Welch test to the case of arbitrarily many samples. Welch’s correction reduces the denominator degrees of freedom in the F-test.
  • Question: I read there are no SS’s in Welch’s ANOVA. Why? *pwr::unknown function but the sjstats package must use it to calculate Cohen’s D!
  • Multiple comparisons (Games-Howell): South not different from Flat!
#1-way Welch's ANOVAs
W.Aovs<- D %>% select(Fv.Fm) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Topography, var.equal= FALSE, data= D)), 
                .id="Variable")#run anova
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract R-sq values from summary.lm(Mod)
R.sq<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-summary.lm(Mod1)
     Out$r.squared
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Welch's ANOVA Fv/Fm vs Toopography - heteroskedastic", digits = c(1, 1, 3, 3, 4, 5, 3), table.attr = "style = \"color: black;\"") %>% kable_styling()
Welch’s ANOVA Fv/Fm vs Toopography - heteroskedastic
Variable num.df den.df statistic p.value Rsq Cohens_F
Fv.Fm 2 54.745 12.098 0 0.1932 0.489

#Multiple comparisons - Games-Howell
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following objects are masked from 'package:plyr':
## 
##     desc, mutate
## The following object is masked from 'package:stats':
## 
##     filter
games_howell_test(data=D, formula= Fv.Fm ~ Topography, conf.level=0.95, detailed=F)

Residual plots

  • Goal: make residual plots of the Welch’s ANOVA model (Fv.Fm ~ Topography) to be fit below to test assumption of normality.
  • Results: QQplot looks pretty normal! Assumption satisfied!
#Fv/Fm ~ Topography
Mod1<- aov(Fv.Fm ~ Topography, data=D)
summary.lm(Mod1)
## 
## Call:
## aov(formula = Fv.Fm ~ Topography, data = D)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.288588 -0.081234 -0.001326  0.095766  0.255541 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.43183    0.02823  15.298   <2e-16 ***
## TopographyF -0.04037    0.03595  -1.123   0.2644    
## TopographyN  0.10676    0.03655   2.921   0.0044 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1354 on 91 degrees of freedom
## Multiple R-squared:  0.1932, Adjusted R-squared:  0.1754 
## F-statistic: 10.89 on 2 and 91 DF,  p-value: 5.739e-05
plot(Mod1)

Classic ANOVA: PhiPSII ~ Topography

  • use stats::oneway.test() with argument var.equal = TRUE, for classic ANOVA.
  • Multiple comparisons result: Phi-PSII differs only between North & Flat!
#1-way Welch's ANOVAs
W.Aovs<- D %>% select(PhiPSII) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Topography, var.equal= TRUE, data= D)), 
                .id="Variable")#run anova
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract R-sq values from summary.lm(Mod)
R.sq<-D %>% select(PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-summary.lm(Mod1)
     Out$r.squared
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Extract Cohen's F from batch
CF<-D %>% select(PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "ANOVA: Phi vs Toopography - equal variances", digits = c(1, 1, 3, 3, 4, 5, 3), table.attr = "style = \"color: black;\"") %>% kable_styling()
ANOVA: Phi vs Toopography - equal variances
Variable num.df den.df statistic p.value Rsq Cohens_F
PhiPSII 2 91 6.483 0.0023 0.1247 0.377

#Multiple comparisons
Mod<- aov(PhiPSII ~ Topography, data = D, var.equal= TRUE)
TukeyHSD(Mod)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = PhiPSII ~ Topography, data = D, var.equal = TRUE)
## 
## $Topography
##            diff         lwr        upr     p adj
## F-S -0.06330905 -0.14529288 0.01867479 0.1624954
## N-S  0.04706138 -0.03629741 0.13042018 0.3741169
## N-F  0.11037043  0.03701940 0.18372146 0.0015646

Habitat ANOVAs

  • Equal variances, so classic anovas:
  • Multiple comparisons (Fv/Fm): Canopy differs from dripline & interspace (interspace = dripline)!
  • Multiple comparisons (Phi): All three habitat types differ!!
#1-way classic ANOVAs - use var.equal = TRUE
W.Aovs<- D %>% select(Fv.Fm, PhiPSII) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Habitat, var.equal= TRUE, data= D)), 
                .id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, denom.df
## Multiple parameters; naming those columns num.df, denom.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Habitat, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Fv.Fm, PhiPSII) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Habitat, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Classic ANOVAs vs Microhabitat type - equal variance", digits = c(3, 3, 3, 3, 4, 5, 4), table.attr = "style = \"color: black;\"") %>% kable_styling()
Classic ANOVAs vs Microhabitat type - equal variance
Variable num.df den.df statistic p.value Cohens_F Rsq
Fv.Fm 2 91 12.433 0 0.523 0.215
PhiPSII 2 91 14.078 0 0.556 0.236

#Multiple comparisons (Fv/Fm)
Mod<- aov(Fv.Fm ~ Habitat, data = D, var.equal= TRUE)
TukeyHSD(Mod)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Fv.Fm ~ Habitat, data = D, var.equal = TRUE)
## 
## $Habitat
##                            diff        lwr         upr     p adj
## Dripline-Canopy     -0.09583631 -0.1715138 -0.02015878 0.0091863
## Interspace-Canopy   -0.17431250 -0.2622704 -0.08635459 0.0000250
## Interspace-Dripline -0.07847619 -0.1746205  0.01766808 0.1321162

#Multiple comparisons (Phi)
Mod<- aov(PhiPSII ~ Habitat, data = D, var.equal= TRUE)
TukeyHSD(Mod)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = PhiPSII ~ Habitat, data = D, var.equal = TRUE)
## 
## $Habitat
##                            diff        lwr         upr     p adj
## Dripline-Canopy     -0.11274107 -0.1813238 -0.04415837 0.0005046
## Interspace-Canopy   -0.15511806 -0.2348298 -0.07540627 0.0000349
## Interspace-Dripline -0.04237698 -0.1295077  0.04475369 0.4807655